This code contains model files and associated material for the 2023 ICES course on Stock Synthesis held in Copenhagen, Denmark.
# install.packages("devtools")
# devtools::install_github("r4ss/r4ss", ref="development")
#
# install.packages("caTools")
# library("caTools")
#
# install.packages("r4ss")
library(r4ss)
library(here)
#remotes::install_github("PIFSCstockassessments/ss3diags")
library(ss3diags)
library(knitr)
library(tidyverse)
dir1<-here("herring-aspm")# version 3.30.21
dir2<-here("Workshops" , "SS3_ICES_Course", "herring-just-2019-as-equilibrium")
dir3<-here("herring-scaa") #
dir4<-here("Workshops" , "SS3_ICES_Course", "herring-scaa-2dar") #
dir5<-here("herring-scaa-2dar-noTV")
dir6<-here("Workshops" , "SS3_ICES_Course", "predators") # v. 3.30.21
Run models
#shell(cmd="ss_win") # run SS para windows
# aqui debo cambiar el setwd() para correr el modelo
#dir5 <- setwd("~/DOCAS/SA_Krill/s5")
# r4ss::run(dirvec = dir3, model = './ss_osx',
# skipfinished = FALSE)
# #system('./ss_osx')
#another way to run with new r4ss
r4ss::run(
dir = dir1,
exe = "ss_osx",
skipfinished = FALSE, # TRUE will skip running if Report.sso present
show_in_console = TRUE # change to true to watch the output go past
)
r4ss::run(
dir = dir1,
exe = "ss_osx",
skipfinished = FALSE, # TRUE will skip running if Report.sso present
show_in_console = TRUE # change to true to watch the output go past
)
r4ss::run(
dir = dir2,
exe = "ss_osx",
skipfinished = FALSE, # TRUE will skip running if Report.sso present
show_in_console = TRUE # change to true to watch the output go past
)
r4ss::run(
dir = dir3,
exe = "ss_osx",
skipfinished = FALSE, # TRUE will skip running if Report.sso present
show_in_console = TRUE # change to true to watch the output go past
)
r4ss::run(
dir = dir4,
exe = "ss_osx",
skipfinished = FALSE, # TRUE will skip running if Report.sso present
show_in_console = TRUE # change to true to watch the output go past
)
r4ss::run(
dir = dir5,
exe = "ss_osx",
skipfinished = FALSE, # TRUE will skip running if Report.sso present
show_in_console = TRUE # change to true to watch the output go past
)
r4ss::run(
dir = dir6,
exe = "ss_osx",
skipfinished = FALSE, # TRUE will skip running if Report.sso present
show_in_console = TRUE # change to true to watch the output go past
)
herring-aspm ModelSaco los outputs en html
Data disponible para este escenario. Espinel es la serie mas consistente del conjunto de datos.
SSplotData(base.model1, subplot = 2)
## Warning: The `subplot` argument of `SSplotData()` is deprecated as of r4ss 1.45.1.
## ℹ Please use the `subplots` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Respecto a los valores y parametros biologicos modelados, los siguientes graficos identifican los estimadores puntuales del recurso
SSplotBiology(base.model,
subplots =1,
labels = c("Length (cm)", "Age (yr)",
"Maturity",
"Mean weight (kg) in last year",
"Spawning output",
"Length (cm, beginning of the year)",
"Natural mortality",
"Female weight (kg)",
"Female length (cm)",
"Fecundity",
"Default fecundity label",
"Year",
"Hermaphroditism transition rate",
"Fraction females by age at equilibrium"),
)
aporte de las cohortes por año para las capturas.
SSplotCohortCatch(base.model1)
AJuste de tallas por flota
SSplotComps(base.model1, subplots = 1)
Otros plots
SSplotDynamicB0(base.model1, uncertainty = F)
## Plotting Dynamic B0
#SSplotSPR(base.model3)
SSplotPars(base.model1)
## Plotting distributions for 4 estimated parameters (deviations not included).
## $mcmc not found in input 'replist', changing input to 'showpost=FALSE'
Salida de las biomasas con las dos flotas
SSplotTimeseries(base.model1, subplot = 1)
SSplotTimeseries(base.model1,
subplot=1,
add = FALSE,
forecastplot = FALSE,
uncertainty = FALSE,
bioscale = 1,
minyr = -Inf,
maxyr = Inf,
plot = TRUE,
print = FALSE,
plotdir = "default",
verbose = TRUE,
btarg = "default",
minbthresh = "default",
xlab = "Year",
labels = NULL,
pwidth = 6.5,
pheight = 5,
punits = "in",
res = 300,
ptsize = 10,
cex.main = 1,
mainTitle = FALSE,
mar = NULL
)
herring-scaa ModelSaco los outputs en html
Data disponible para este escenario. Espinel es la serie mas consistente del conjunto de datos.
SSplotData(base.model3, subplot = 2)
SSplotRecdevs(base.model3)
Respecto a los valores y parametros biologicos modelados, los siguientes graficos identifican los estimadores puntuales del recurso
SSplotBiology(base.model3, subplots =1,
labels = c("Length (cm)",
"Age (yr)",
"Maturity",
"Mean weight (kg) in last year",
"Spawning output",
"Length (cm, beginning of the year)",
"Natural mortality",
"Female weight (kg)",
"Female length (cm)",
"Fecundity",
"Default fecundity label",
"Year",
"Hermaphroditism transition rate",
"Fraction females by age at equilibrium"),
)
aporte de las cohortes por año para las capturas.
SSplotCohortCatch(base.model3)
AJuste de tallas por flota
SSplotComps(base.model3, subplots = 1)
Otros plots
SSplotDynamicB0(base.model3, uncertainty = F)
## Plotting Dynamic B0
#SSplotSPR(base.model3)
SSplotPars(base.model3)
## Excluding 57 deviation parameters because input 'showdev' = FALSE
## Plotting distributions for 14 estimated parameters (deviations not included).
## $mcmc not found in input 'replist', changing input to 'showpost=FALSE'
Salida de las biomasas con las dos flotas
SSplotTimeseries(base.model3, subplot = 10)
SSplotTimeseries(base.model3,
subplot=1,
add = FALSE,
forecastplot = FALSE,
uncertainty = FALSE,
bioscale = 1,
minyr = -Inf,
maxyr = Inf,
plot = TRUE,
print = FALSE,
plotdir = "default",
verbose = TRUE,
btarg = "default",
minbthresh = "default",
xlab = "Year",
labels = NULL,
pwidth = 6.5,
pheight = 5,
punits = "in",
res = 300,
ptsize = 10,
cex.main = 1,
mainTitle = FALSE,
mar = NULL
)
herring-scaa-2dar-noTV ModelSaco los outputs en html
Data disponible para este escenario. Espinel es la serie mas consistente del conjunto de datos.
SSplotData(base.model5, subplot = 2)
Respecto a los valores y parametros biologicos modelados, los siguientes graficos identifican los estimadores puntuales del recurso
SSplotBiology(base.model5, subplots =1,
labels = c("Length (cm)",
"Age (yr)",
"Maturity",
"Mean weight (kg) in last year",
"Spawning output",
"Length (cm, beginning of the year)",
"Natural mortality",
"Female weight (kg)",
"Female length (cm)",
"Fecundity",
"Default fecundity label",
"Year",
"Hermaphroditism transition rate",
"Fraction females by age at equilibrium"),
)
## Note: this model uses the empirical weight-at-age input.
## Plots of many quantities related to growth are skipped.
## Warning in max(biology[["Len_mean"]]): no non-missing arguments to max;
## returning -Inf
aporte de las cohortes por año para las capturas.
SSplotCohortCatch(base.model5)
AJuste de tallas por flota
SSplotComps(base.model5, subplots = 1)
Otros plots
SSplotDynamicB0(base.model5, uncertainty = F)
## Plotting Dynamic B0
#SSplotSPR(base.model3)
SSplotPars(base.model5)
## Excluding 69 deviation parameters because input 'showdev' = FALSE
## Plotting distributions for 3 estimated parameters (deviations not included).
## $mcmc not found in input 'replist', changing input to 'showpost=FALSE'
Salida de las biomasas con las dos flotas
SSplotTimeseries(base.model5, subplot = 10)
SSplotTimeseries(base.model5,
subplot=1,
add = FALSE,
forecastplot = FALSE,
uncertainty = FALSE,
bioscale = 1,
minyr = -Inf,
maxyr = Inf,
plot = TRUE,
print = FALSE,
plotdir = "default",
verbose = TRUE,
btarg = "default",
minbthresh = "default",
xlab = "Year",
labels = NULL,
pwidth = 6.5,
pheight = 5,
punits = "in",
res = 300,
ptsize = 10,
cex.main = 1,
mainTitle = FALSE,
mar = NULL
)
Los desembarques utilizados para cada una de las flotas que generan remoción en el recurso reineta, a saber; espinel, enmalle e industrial (Figura 6).
Desembarques de merluza común por flota
Random Walk (RW) refers to a mathematical model that describes a stochastic process in which a variable changes randomly over time, without a clear trend or pattern.
Specifically, a random walk can be used as a Bayesian estimation technique to infer the posterior distribution of an unknown parameter. In this approach, it is assumed that the prior distribution of the parameter is a normal distribution with mean zero and a known variance, and that the parameter value at each time step follows a random walk process. Based on the observed data and the prior distribution, the posterior distribution of the parameter can be calculated using Bayesian inference.
Random walk is a useful tool for parameter estimation in dynamic models, as it allows for modeling uncertainty and variability in the parameter’s evolution over time. However, it is important to note that the random walk assumes that the changes in the parameter are random and without a clear trend, which may not be appropriate in all cases.
#########################################################################
#COMPARACION DE MODELOS con distinto desembarque
#########################################################################
#PLOT labels, name of each model run
legend.labels <- c('herring-aspm',
'herring-just-2019-as-equilibrium')
#read in all model runs
#note if cover=T you need a hessian; if covar=F you do not need a hessian
biglist <- SSgetoutput(keyvec = NULL,
dirvec = c(dir1,dir2,dir3),
getcovar = F)
#create summary of model runs from list above
summaryoutput <- SSsummarize(biglist)
SSplotComparisons(summaryoutput, subplots = 1:20, plot = TRUE,
print = T, endyrvec = "default", indexfleets = NULL,
indexUncertainty = FALSE,
indexQlabel = TRUE,
indexQdigits = 4,
indexSEvec = "default",
indexPlotEach = FALSE,
labels = c("Year", "Spawning biomass (t)",
"Relative spawning biomass",
"Age-0 recruits (1,000s)",
"Recruitment deviations",
"Index",
"Log index",
"1 - SPR",
"Density",
"Management target",
"Minimum stock size threshold",
"Spawning output",
"Harvest rate"), col = NULL,
shadecol = NULL, pch = NULL,
lty = 1, lwd = 2, spacepoints = 10, staggerpoints = 1,
initpoint = 0, tickEndYr = TRUE, shadeForecast = TRUE,
xlim = "default", ylimAdj = 1, xaxs = "r", yaxs = "r",
type = "o", uncertainty = TRUE,
shadealpha = 0.1, legend = TRUE,
legendlabels = "default", legendloc = "topright",
legendorder = "default", legendncol = 1, sprtarg = NULL,
btarg = NULL, minbthresh = NULL,
pwidth = 6.5, pheight = 5,
punits = "in", res = 300, ptsize = 10,
plotdir = "~/DOCAS/Workshops/SS3_ICES_Course",
filenameprefix = "",
densitynames = c("SSB_Virgin", "R0"),
densityxlabs = "default", densityscalex = 1,
densityscaley = 1, densityadjust = 1,
densitysymbols = TRUE,
densitytails = TRUE,
densitymiddle = FALSE,
densitylwd = 1,
fix0 = TRUE, new = TRUE,
add = FALSE,
par = list(mar = c(5, 4, 1, 1) + 0.1),
verbose = TRUE, mcmcVec = FALSE,
show_equilibrium = TRUE)
SStableComparisons(summaryoutput,
likenames = c("TOTAL",
"Survey",
"Length_comp",
"Age_comp",
"priors",
"Size_at_age"),
names = c("Recr_Virgin",
"R0", "steep",
"NatM",
"L_at_Amax",
"VonBert_K",
"SSB_Virg",
"Bratio_2017", "SPRratio_2016"),
digits = NULL, modelnames = "default", csv = FALSE,
csvdir ="~/DOCAS/Workshops/SS3_ICES_Course",
csvfile = "parameter_comparison_table.csv", verbose = TRUE,
mcmc = FALSE)
#do retrospective model runs
SS_doRetro('~/DOCAS/Workshops/SS3_ICES_Course',
'', newsubdir = "retrospectives",
subdirstart = "retro", years = 0:-5,
overwrite = TRUE, exefile = "ss",
extras = "-nox -nohess",
intern = FALSE, CallType = "system",
RemoveBlocks = FALSE)
retroModels <- SSgetoutput(dirvec=file.path('~/DOCAS/Workshops/SS3_ICES_Course',
"retrospectives",
paste("retro",0:-5,sep="")))
retroSummary <- SSsummarize(retroModels)
endyrvec <- retroSummary$endyrs + 0:-5
SSplotComparisons(retroSummary,
endyrvec=endyrvec,
legendlabels=paste("Data",0:-5,"years"),
plotdir='~/DOCAS/Workshops/SS3_ICES_Course',
plot=TRUE,print=T)
TableCompare <- SStableComparisons(retroSummary,
likenames=like,
names=names,
modelnames=c('B','-1','-2','-3','-4','-5'),
csv=TRUE,
csvfile="RetroRuns.csv",verbose=FALSE)